﻿class correlador
{
private:
public:
	//Longitud de los correladores
	int ncorr;
	//Tiempo
	double *tiempo;
	//Desplazamiento cuadratico medio
	double *msd;
	//Correladores del calculo de la viscosidad
	double *A12,*A13,*A23;
	//Conteo para la normalizacion
	int *jcorr;
	//
	double *dn;
	//Tiempos de relajacion
	double *trelaj;
	//Almacenamiento posiciones particulas para msd
	double **x0c, **y0c,**z0c;
	//Valor para el tiempo de relajacion del primer canal
	double trelaj1;
	//Caja
	caja *box;
	//Paso de tiempo
	double dt, dth;
	//Determina la viscosidad aprox = (beta/6V)*lim t-> inf (A*A/t) 	
	double sumtau;
	//Solo componentes xz e yz del tensor
	double sumtau2;
	//Nombre del fichero
	parametros paramet;
	//Anchura de los campos de salida
	int ancho;

	correlador(caja *bx, parametros param)
	{
		ancho=20;
		paramet = param;
		ncorr=param.ncorr;
		box=bx;
		trelaj1=param.trelaj1;
		dt=param.Dt;
		dth=0.5*dt;
		//Contruimos los arrays
		tiempo = new double[ncorr];
		msd = new double[ncorr];
		A12 = new double[ncorr];
		A13 = new double[ncorr];
		A23 = new double[ncorr];
		jcorr = new int[ncorr];
		dn = new double[ncorr];
		trelaj = new double[ncorr];
		x0c = new double*[ncorr];
		y0c = new double*[ncorr];
		z0c = new double*[ncorr];
		for(int i=0;i<ncorr;i++)
		{
			x0c[i]=new double[box->npart];
			y0c[i]=new double[box->npart];
			z0c[i]=new double[box->npart];
		}

		inicializa();
	}

	~correlador()
	{
		//Eliminamos los arrays
		delete []tiempo;
		delete []msd;
		delete []A12;
		delete []A13;
		delete []A23;
		delete []jcorr;
		delete []dn;
		delete []trelaj;

		for(int i=0;i<box->npart;i++)
		{
			delete x0c[i];
			delete y0c[i];
			delete z0c[i];
		}
		delete []x0c;
		delete []y0c;
		delete []z0c;
	}

	void inicializa(void);
	void calcCorrelador(double t);
	//Inicializacion del archivo de correlacion
	void initFileCorr(void);

};

void correlador::inicializa(void)
{
	for(int icorr=0; icorr<ncorr;icorr++)
	{
		dn[icorr]=1.0;
		jcorr[icorr]=0;
		A12[icorr]=0.0;
		A13[icorr]=0.0;
		A23[icorr]=0.0;
		if (icorr==0)
		{
			//Valor del primer canal del tiempo de relajacion
			trelaj[icorr]=trelaj1;
		}
		else
		{
			//Se les da a todos los correladores un tiempo de relajacion enorme excepto al primero
			trelaj[icorr]=1.0e8;
		}
	}
}

void correlador::initFileCorr(void)
{
	ofstream MiArchivo (paramet.corrFile,ios::out);
	MiArchivo.precision(4);
	if (MiArchivo.is_open())
	{
		MiArchivo.setf(ios::fixed);
		MiArchivo<<"Correladores"<<endl;
		MiArchivo.fill('*');
		MiArchivo << setw(8*ancho) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<setw(ancho)<<"icorr"<<setw(ancho)<<"t-trelaj"<<setw(ancho)<<"dcmijx"<<setw(ancho)<<"dcmijy"<<setw(ancho)<<"dcmijz"
				<<setw(ancho)<<"sumtau"<<setw(ancho)<<"sumtau2"<<endl;
		MiArchivo.fill('-');
		MiArchivo << setw(8*ancho) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo.close();
	}

}

void correlador::calcCorrelador(double t)
{

	double dd[3];
	double dcmij[3];
	double dcm;
	double dummy;
	int idn;

	ofstream MiArchivo (paramet.corrFile,ios::out|ios::app);
	MiArchivo.precision(4);

	//Barre todos los correladores que ya han empezado a medir
	for(int icorr=0;icorr<ncorr;icorr++)
	{
		//dn(icorr) es el salto de tiempo del correlador icorr con el que escribe 
		//Conforme el correlador tal va midiendo, lo hace con menos frecuencia cada 
		//vez. Es decir, al principio mide mas seguido. Conforme pasa el tiempo, una 
		//medida y otra van espaciandose cada vez mas. Se alarga el paso.
		idn=int(dn[icorr]) ;
		//Si el tiempo por el que va el programa supera al tiempo de relajaci�n de tal
		//correlador entonces ese correlador empieza a medir y se va contando las 
		//veces que ese correlador aparece usando jcorr
		if(t>=trelaj[icorr]-dth)
		{

			jcorr[icorr]++;

			//Integrales 1,2; 1,3 y 2,3 por el metodo del rectengulo. alfa<beta 
			//Se va integrando para calcular finalmente la viscosidad

			A12[icorr]+=(box->spxy)*dt;
			A13[icorr]+=(box->spxz)*dt;
			A23[icorr]+=(box->spyz)*dt;

			//Si esta condicion se cumple quiere decir que estamos en 
			//la primera medida del correlador tal

			if (abs(t-trelaj[icorr])<dth)
			{
				//Se asignan las posiciones del instante actual a las posiciones 
				//de las que parte el correlador tal como sus posiciones iniciales. 
				//Se hace para todas las particulas

				for(int k=0;k<(box->npart);k++)
				{
					x0c[icorr][k]=box->part[k].pos[0];
					y0c[icorr][k]=box->part[k].pos[1];
					z0c[icorr][k]=box->part[k].pos[2];
				}
			}

			//Si el numero de veces que ha contado el correlador tal es mayor que 
			//la parte entera de la frecuencia con la que esta midiendo, escribe el 
			//tiempo por el que va menos el tiempo de relajacion del correlador tal
			//y el correlador que es. En este momento es cuando se calcula el DCM 
			//para cada eje y se escriben las pos y vel de las 1000 particulas en 
			//el archivo "pos.dat"

			if (jcorr[icorr]>=idn)
			{
				//write(12) t-trelaj(icorr),icorr
				dcm=0.;
				for(int k=0;k<3;k++)
				{
					dcmij[k]=0.0;
				}

				for(int j=0;j<box->npart;j++)
				{
					//Escribe posiciones y velocidades
					//write (12) x0(j),y0(j),z0(j),vx0(j),vy0(j),vz0(j) 

					//Correccion de la distancia entre las particulas i y j

					dd[0]=box->part[j].pos[0]-x0c[icorr][j];                          
					dd[0]-=(box->lbox[0])*anint(dd[0]/(box->lbox[0])); //Entero + proximo manteniendo signo 
					dd[1]=box->part[j].pos[1]-y0c[icorr][j];
					dd[1]-=(box->lbox[1])*anint(dd[1]/(box->lbox[1]));
					dd[2]=box->part[j].pos[2]-z0c[icorr][j];
					dd[2]-=(box->lbox[2])*anint(dd[2]/(box->lbox[2]));

					//Calculo de los desplazamientos cuadraticos medios
					for(int k=0;k<3;k++)
					{
						dummy = dd[k]*dd[k];
						dcmij[k]+=dummy;
						dcm+=dummy;
					}
				}
				for(int k=0;k<3;k++)
				{
					dcmij[k]/=double(box->npart);
				}
				dcm/=double(box->npart);

				//Determina la viscosidad aprox = (beta/6V)*lim t-> inf (A*A/t) 
				dummy = (A13[icorr]*A13[icorr]+A23[icorr]*A23[icorr]);
				sumtau=(A12[icorr]*A12[icorr]+dummy)/3.;
				//Solo componentes xz e yz del tensor
				sumtau2=(dummy)/2.; 

				if((t-trelaj[icorr])>1.0e-5)
				{
					if (MiArchivo.is_open())
					{
						//escribe cuando ha pasado el tiempo de relajacion del correlador icorr
						MiArchivo<<setw(ancho)<<icorr<<setw(ancho)<<t-trelaj[icorr];
						for(int k=0;k<3;k++)
						{
							MiArchivo<<setw(ancho)<<dcmij[k];
						}
						MiArchivo<<setw(ancho)<<sumtau/(t-trelaj[icorr])<<setw(ancho)<<sumtau2/(t-trelaj[icorr])<<endl;
					}
				}

				if (icorr<ncorr)
				{
					//si el correlador por el que va mide un DCM mayor de 0.5 y ademas 
					//el siguiente correlador aun no se ha iniciado, iguala el tiempo 
					//por el que va el programa con el tiempo de relajacion del 
					//correlador siguiente para que en la proxima ronda empiece 
					//a medir este ultimo

					if ((dcm>0.5)&&(trelaj[icorr+1]==1e8))
					{
						trelaj[icorr+1]=t;
					}
				}

				//Va incrementando el paso para escribir
				dn[icorr]*=1.25 ;
			}

			//Si se han cubierto todos los correladores y se ha superado un DCM de 25, 
			//entonces se acaba la ejecucion
			if ((icorr==ncorr) && (dcm>25.)) 
			{
				break; 
			}
		}

	}
	if (MiArchivo.is_open())
	{
		MiArchivo.close();
	}
}